/*===================================================================
analysis.do
	*Author(s):	David Phillips, Sean McConville, Charlie Law
	*Purpose:	Uses processed data to make final tables and figures
===================================================================*/

	*do "$mydir/Scheme_LEO_v2"

	cd "$outputdir"

	
***
***Main effects (Tables 2, A.6, A.7, A.8)
***
	
	*looping across samples and specifications
	foreach samp in early late controls veryearly clustered {
		clear
		use "$datadir/regdata.dta"
		gen controls = early
		gen clustered = early
		keep if `samp' == 1
		
		if "`samp'" ~= "controls" {
			global controlvars 
		}		
		
		if "`samp'" == "controls" {
			global controlvars $spdatlist1  $spdatlist2 $lagoutcomelist new_add_pre_mis
			gen new_add_pre_mis = (new_addpre_12 == .)
			replace new_addpre_12 = 0 if new_addpre_12 == .
		}

		if "`samp'" ~= "clustered" {
			global clustvar HMISID
		}

		if "`samp'" == "clustered" {
			global clustvar strata
		}
		
				
		
		***  Outcome table in balance form
			
			gen dummy = .
			gen const=1 

			
			forvalues x = 1(1)12 {															
				label var dh_`x'm "D:H Homelessness Prevention (`x' mos)"
				label var hp_`x'mo "Any Homelessness Prevention (`x' mos)"
				label var non_hp_`x'm "Homeless (`x' mos)"
				label var non_hp3_`x'm "Homeless (3 to `x' mos)"
				label var new_add_`x' "Address Change (`x' mos)"
				label var new_add_notscc_`x' "Move Out of County (`x' mos)"
				label var shelter_`x'm "--Shelter (`x' mos)"
				label var perm_`x'm "--Permanent Supportive Housing (`x' mos)"
				label var tran_`x'm "--Transitional Housing (`x' mos)"
				label var rrh_`x'm "--Rapid Rehousing (`x' mos)"
				label var out_`x'm "--Outreach (`x' mos)"
				label var coord_`x'm "--Coordinated Entry (`x' mos)"
				label var big_`x'm "--Other Homeless Services (`x' mos)"
				label var paymtany_tot_`x' "Any Payment (`x' mos)"
				label var paymt_tot_`x' "Total Payments (`x' mos)"
				label var paymt_rent_`x' "--Rent (`x' mos)"
				label var paymt_sd_`x' "--Security Deposit (`x' mos)"
				label var paymt_util_`x' "--Utilities (`x' mos)"
				label var paymt_other_`x' "--Other (`x' mos)"

			}
	
			forvalues x = 1(1)12 {															
				label var dh_`x'm_days "D:H Homelessness Prevention Days (`x' mos)"
				label var hp_`x'mo_days "Homelessness Prevention Days (`x' mos)"
				label var non_hp_`x'm_days "Homeless Days (`x' mos)"
				label var new_add_`x' "Address Change (`x' mos)"
				label var shelter_`x'm_days "Shelter Days (`x' mos)"
				label var perm_`x'm_days "Permanent Supportive Housing Days (`x' mos)"
				label var tran_`x'm_days "Transitional Housing Days (`x' mos)"
				label var rrh_`x'm_days "Rapid Rehousing Days (`x' mos)"
				label var out_`x'm_days "Outreach Days (`x' mos)"
				label var coord_`x'm_days "Coordinated Entry Days (`x' mos)"
				label var big_`x'm_days "Subsidized Housing (`x' mos)"
			}
	
			
			reg dh_3m const
			eststo _filler
			
			***First Stage
			
				esttab _filler _filler _filler _filler using "Tables\hp_outcomesnoadd_`samp'.tex", replace ///
				keep(const) drop(const) ///
				mtitles("Control" "Treatment" "Dif (OLS)" "Dif (IV)") ///
				nonotes noobs nogaps fragment 
				
				foreach var of varlist dh_3m hp_3mo paymtany_tot_3 paymt_tot_3 paymt_rent_3 paymt_sd_3 paymt_other_3 {

					replace dummy=const
					
					qui reg `var' dummy if treatment_sum == 0   , nocons robust
						eststo _reg1
						
					qui reg `var' dummy if treatment_sum == 1  , nocons robust
						eststo _reg2		
					
					replace dummy=treatment_sum
					
					qui reg `var' dummy i.strata $controlvars , cluster($clustvar)
						eststo _reg3
						
					replace dummy = hp_3mo
					qui xi: ivreg `var'  (dummy	 = treatment_sum)  i.strata $controlvars , cluster($clustvar)
						eststo _reg4
							
							
					qui su `var'
					scalar varmean = r(mean)
					forvalues x = -10(1)10 {
						if varmean > 10^`x' {
							scalar digits = `x' 
						}
					}
					scalar digits = max(-digits + 1,0)
	
					esttab _reg1 _reg2  _reg3  using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep(dummy) varlabels(dummy "`: var label `var''") ///
					cells("b(fmt(%9.`=`=digits''f) pattern(1 1 0 0)) b(fmt(%9.`=`=digits''f) star pattern(0 0 1))") ///
					starlevel(* 0.10 ** 0.05 *** 0.01) ///
					collabels(none) nomtitles ///
					nodepvars nonum noobs nonotes nolines fragment
					
					esttab  _reg1 _reg2  _reg3 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep(dummy) varlabels(dummy " ") ///
					cells(se(fmt(%9.`=`=digits''f) par pattern(0 0 1)) ) ///
					starlevel(* 0.10 ** 0.05 *** 0.01) ///
					collabels(none) nomtitles ///
					nodepvars nonum noobs nonotes b(a4) fragment nolines					
					
				}
				
			
				esttab  _reg1 _reg2  _reg3 _reg4 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
					nonotes fragment nolines nodepvars nonotes nogaps
					
 				esttab  _reg1 _reg2  _reg3 _reg4 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
 					keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
 					nonotes fragment noobs
					
								
			***HMIS Outcomes
				
				global hmisvars non_hp_6m  shelter_6m out_6m big_6m  non_hp_6m_days   
				if  "`samp'" == "early" | "`samp'" == "controls" {
					global hmisvars non_hp_6m  shelter_6m out_6m big_6m  non_hp_12m non_hp3_12m non_hp_6m_days  non_hp_12m_days  
				}
								
				foreach var of varlist $hmisvars  {

					replace dummy=const
					
					qui reg `var' dummy if treatment_sum == 0   , nocons robust
						eststo _reg1
						
					qui reg `var' dummy if treatment_sum == 1  , nocons robust
						eststo _reg2		
					
					replace dummy=treatment_sum
					
					qui reg `var' dummy i.strata  $controlvars , cluster($clustvar)
						eststo _reg3
						
					replace dummy = hp_3mo
					qui xi: ivreg `var'  (dummy	 = treatment_sum)  i.strata $controlvars , cluster($clustvar)
						eststo _reg4

					qui su `var'
					scalar varmean = r(mean)
					forvalues x = -10(1)10 {
						if varmean > 10^`x' {
							scalar digits = `x' 
						}
					}
					scalar digits = max(-digits + 1,0)
		
					esttab _reg1 _reg2  _reg3 _reg4 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep(dummy) varlabels(dummy "`: var label `var''") ///
					cells("b(fmt(%9.`=`=digits''f) pattern(1 1 0 0)) b(fmt(%9.`=`=digits''f) star pattern(0 0 1 1))") ///
					starlevel(* 0.10 ** 0.05 *** 0.01) ///
					collabels(none) nomtitles ///
					nodepvars nonum noobs nonotes nolines fragment
					
					esttab  _reg1 _reg2  _reg3 _reg4 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep(dummy) varlabels(dummy " ") ///
					cells(se(fmt(%9.`=`=digits''f) par pattern(0 0 1 1)) ) ///
					starlevel(* 0.10 ** 0.05 *** 0.01) ///
					collabels(none) nomtitles ///
					nodepvars nonum noobs nonotes b(a4) fragment nolines		
					
				}
				
					esttab  _reg1 _reg2  _reg3 _reg4 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
					nonotes fragment nolines nogaps
					
					esttab  _reg1 _reg2  _reg3 _reg4 using "Tables\hp_outcomesnoadd_`samp'.tex", append ///
					keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
					nonotes fragment noobs
					
		
		}



		*test that address change and HMIS effects are different
		
		clear
		use "$datadir/regdata.dta"
		keep if early == 1
		
		reg non_hp_6m treatment_sum i.strata
		egen num = total(1), by(strata)
		*avoiding singular covariance matrix for strata with small sample sizes; doesn't change estimates barely at all
		replace strata = 999 if num < 3
		*
		reg non_hp_6m treatment_sum i.strata
		eststo _reg1
		reg new_add_6  treatment_sum i.strata
		eststo _reg2
		suest _reg1 _reg2, robust
		nlcom _b[_reg1_mean:treatment_sum] - _b[_reg2_mean:treatment_sum] 
			
		tab non_hp_6m new_add_6, col
		





***
*** More detailed payments table (Table A.5) 
***

	foreach samp in  early  {
		clear
		use "$datadir/regdata.dta"
		gen controls = early
		keep if `samp' == 1
			
		foreach var in tot rent sd util trans other_hous other {
			gen paymtcop_`var'_3 = paymt_`var'_3
			replace paymtcop_`var'_3 = . if paymt_`var'_3 == 0
		}

		gen dummy = .
		gen const = 1

			forvalues x = 1(1)12 {																	
				label var paymtany_tot_`x' "Any Payment (`x' mos)"
				label var paymtany_rent_`x' "--Any Rent (`x' mos)"
				label var paymtany_sd_`x' "--Any Security Deposit (`x' mos)"
				label var paymtany_other_`x' "--Any Other (`x' mos)"
				label var paymtany_other_hous_`x' "----Any Other Housing (`x' mos)"
				label var paymtany_util_`x' "----Any Utilities (`x' mos)"
				label var paymtany_trans_`x' "----Any Transportation (`x' mos)"
				label var paymt_tot_`x' "Total Payments (`x' mos)"
				label var paymt_rent_`x' "--Rent (`x' mos)"
				label var paymt_sd_`x' "--Security Deposit (`x' mos)"
				label var paymt_util_`x' "--Utilities (`x' mos)"
				label var paymt_other_`x' "--Other (`x' mos)"
				label var paymtcop_tot_3 "Total Payment, Conditional on Positive (3 mos)"
				label var paymtcop_rent_3 "--Rent Payment, Conditional on Positive (3 mos)"
				label var paymtcop_sd_3 "--SD Payment, Conditional on Positive (3 mos)"
				label var paymtcop_other_3 "--Other Payment, Conditional on Positive (3 mos)"

			}

			***First Stage
			
				
				esttab _filler _filler _filler using "Tables\hp_paymt_`samp'.tex", replace ///
				keep(const) drop(const) ///
				mtitles("Control" "Treatment" "Dif (OLS)" ) ///
				nonotes noobs nogaps fragment 
				

				foreach var of varlist 	paymt_tot_3 paymt_rent_3 paymt_sd_3 paymt_other_3 ///
										paymtcop_tot_3 paymtcop_rent_3 paymtcop_sd_3 paymtcop_other_3 ///	
										paymtany_tot_3 paymtany_rent_3 paymtany_sd_3 paymtany_other_3 paymtany_util_3 paymtany_other_hous_3 paymtany_trans_3 ///
										{

					replace dummy=const
					
					qui reg `var' dummy if treatment_sum == 0   , nocons robust
						eststo _reg1
						
					qui reg `var' dummy if treatment_sum == 1  , nocons robust
						eststo _reg2		
					
					replace dummy=treatment_sum
					
					qui reg `var' dummy i.strata $controlvars , robust
						eststo _reg3
						
					replace dummy = hp_3mo
					qui xi: ivreg `var'  (dummy	 = treatment_sum)  i.strata $controlvars , robust
						eststo _reg4
							
							
					qui su `var'
					scalar varmean = r(mean)
					forvalues x = -10(1)10 {
						if varmean > 10^`x' {
							scalar digits = `x' 
						}
					}
					scalar digits = max(-digits + 1,0)
						
					esttab _reg1 _reg2  _reg3  using "Tables\hp_paymt_`samp'.tex", append ///
					keep(dummy) varlabels(dummy "`: var label `var''") ///
					cells("b(fmt(%9.`=`=digits''f) pattern(1 1 0 0)) b(fmt(%9.`=`=digits''f) star pattern(0 0 1))") ///
					starlevel(* 0.10 ** 0.05 *** 0.01) ///
					collabels(none) nomtitles ///
					nodepvars nonum noobs nonotes nolines fragment
					
					esttab  _reg1 _reg2  _reg3 using "Tables\hp_paymt_`samp'.tex", append ///
					keep(dummy) varlabels(dummy " ") ///
					cells(se(fmt(%9.`=`=digits''f) par pattern(0 0 1)) ) ///
					starlevel(* 0.10 ** 0.05 *** 0.01) ///
					collabels(none) nomtitles ///
					nodepvars nonum noobs nonotes b(a4) fragment nolines					
					
					
				}
				
				esttab  _reg1 _reg2  _reg3  using "Tables\hp_paymt_`samp'.tex", append ///
					keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
					nonotes fragment nolines nodepvars nonotes nogaps
					
 				esttab  _reg1 _reg2  _reg3  using "Tables\hp_paymt_`samp'.tex", append ///
 					keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
 					nonotes fragment noobs
					
}


***
*** Balance Table (Tables A.1, A.2, A.3)
***

	gen infmatch_early = infmatch * early


	*looping across samples
	foreach samp of varlist main early late infmatch_early {
		clear
		use "$datadir/regdata.dta"
		gen infmatch_early = infmatch * early
		keep if `samp' == 1

		gen dummy = .
		gen const=1 
		

	*making the table
		esttab _filler _filler _filler  using "Tables\hp_balance_`samp'.tex", replace ///
			keep(const) drop(const) ///
			mtitles("Control" "Treat" "Adj Diff") ///
			nonotes noobs  nogaps fragment

			foreach var of varlist $spdatlist1  $spdatlist2 $lagoutcomelist yhat  { 

				replace dummy=const
				
				qui reg `var' dummy if treatment_sum == 0   , nocons robust
					eststo _reg1
				
				qui reg `var' dummy if treatment_sum == 1 , nocons robust
					eststo _reg2		
							
				replace dummy=treatment_sum
				
				qui reg `var' dummy i.strata  , robust
					eststo _reg5
				
				qui su `var'
				scalar varmean = r(mean)
				forvalues x = -10(1)10 {
					if varmean > 10^`x' {
						scalar digits = `x' 
					}
				}
				scalar digits = max(-digits + 1,0)
				
				if "`var'" == "assmtscorez" {
					scalar digits = digits + 1
				}
				
		
				esttab _reg1 _reg2  _reg5  using "Tables\hp_balance_`samp'.tex", append ///
				keep(dummy) varlabels(dummy "`: var label `var''") ///
				cells("b(fmt(%9.`=`=digits''f) pattern(1 1 0)) b(fmt(%9.`=`=digits''f) star pattern(0 0 1))") ///
				starlevel(* 0.10 ** 0.05 *** 0.01) ///
				collabels(none) nomtitles ///
				nodepvars nonum noobs nonotes nolines fragment
				
			}
			
		
				esttab _reg1 _reg2  _reg5 using "Tables\hp_balance_`samp'.tex", append ///
				keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
				nonotes fragment

	}
	
	

	
***
***Context table (Table 1)
***
	clear
	use "$datadir/regdata.dta"

	***table
	gen sample1=0
	replace sample1=1 
	
	gen sample2=0
	replace sample2=1 if ((first_assessment_dt>=td(01jul2019) & first_assessment_dt<=td(31dec2020)) ///
	& sample1==1)
	
	gen sample3=0
	replace sample3=1 if ((assmtscorez<=13 & assmtscorez>=8) & sample2==1)
	
	gen sample4=0
	replace sample4=1 if (( assmtscorez>= 14) & sample2==1)
			
	gen sample5=0
	replace sample5=1 if treatment_sum!=. & eligible_other == 0
	
	gen sample6 = 0
	replace sample6=1 if treatment_sum!=. & eligible_other == 0	& early == 1
	
	gen const = 1
	gen dummy = .
		reg age const
		eststo _filler

	
	*table
	esttab _filler _filler _filler _filler  _filler _filler  _filler  _filler  using "Tables\hp_context.tex", replace ///
		keep(const) drop(const) ///
		mtitles("All" "Always Treat" "RCT Eligible" "In RCT" "In RCT" "In RCT" "In RCT" "In RCT" "In RCT"  ) ///
		nonotes noobs nogaps fragment nolines
		
	esttab _filler _filler _filler _filler _filler _filler _filler _filler using "Tables\hp_context.tex", append ///
		keep(const) drop(const) ///
		mtitles(""  "Score $>$ 13" "8 $\leq$ S $\leq$ 13" "" "Early" "Treat" "Control" "Dif.") ///
		nonotes noobs nonum nogaps fragment

		*foreach var of varlist homeless homeless_times  {
		
		foreach var of varlist $spdatlist $lagoutcomelist yhat {

			replace dummy=const
			
			reg `var' dummy if sample1==1, nocons robust
				eststo _reg1
			
			reg `var' dummy if sample2==1, nocons robust
				eststo _reg2
			
			reg `var' dummy if sample3==1, nocons robust
				eststo _reg3
			
			reg `var' dummy if sample4==1, nocons robust
				eststo _reg4

			reg `var' dummy if sample5==1, nocons robust
				eststo _reg5
				
			reg `var' dummy if sample6==1, nocons robust
				eststo _reg6	

			reg `var' dummy if sample6==1 & treatment_sum == 1, nocons robust
				eststo _reg7
				
			reg `var' dummy if sample6==1 & treatment_sum == 0, nocons robust
				eststo _reg8
				
			replace dummy=treatment_sum
				
			qui reg `var' dummy i.strata  if sample6==1 , robust
				eststo _reg9
					
				
			qui su `var'
			scalar varmean = r(mean)
			forvalues x = -10(1)10 {
				if varmean > 10^`x' {
					scalar digits = `x' 
				}
			}
			scalar digits = max(-digits + 1,0)
				
			esttab _reg2 _reg4 _reg3 _reg5 _reg6 _reg7 _reg8 _reg9 using "Tables\hp_context.tex", append ///
			keep(dummy) varlabels(dummy "`: var label `var''") ///
			cells("b(fmt(%9.`=`=digits''f) pattern(1 1 1 1 1 1 1 0)) b(fmt(%9.`=`=digits''f) star pattern(0 0 0 0 0 0 0 1))") ///
			starlevel(* 0.10 ** 0.05 *** 0.01) ///
			collabels(none) nomtitles ///
			nodepvars nonum noobs nonotes nolines fragment
			
		}

		gen new_addpre_mis = (new_addpre_12 == .)
		replace new_addpre_12 = -99 if new_addpre_12 == .

		reg treatment_sum $spdatlist $lagoutcomelist new_addpre_mis i.strata  if sample6==1 , robust 
			eststo _reg9
			test $spdatlist $lagoutcomelist new_addpre_mis 
			estadd scalar p_value = r(p)

		esttab _reg2 _reg4 _reg3 _reg5 _reg6  _reg7 _reg8 _reg9  using "Tables\hp_context.tex", append ///
			keep (dummy) drop(dummy) nomtitles nonum b(a2) ///
			stats(p_value N, label("Joint test of balance (p-value)" "N") ) //////
			nonotes fragment

		replace new_addpre_12 = . if new_addpre_12 == -99	
			
			

			
			
***
*** Compliers comparison table (Table A.4)
***

	clear
	use "$datadir/regdata.dta"
	keep if early == 1
		
 	global randvar treatment_sum
 	global gottreat hp_3mo
		
	capture file close outfile
	file open outfile using "Tables\ccm.tex", write replace
	file write outfile " & (1)  & (2)  & (3)  & (4) & (5) \\ " _n
	file write outfile " & All  & Compliers  & Always  & Never & (2) - (1) \\ " _n
	file write outfile " & Mean & Mean & Mean & Mean & Dif. \\ \hline " _n
		
		foreach var of varlist $spdatlist $lagoutcomelist yhat  {
		*foreach var of varlist currentbalancz  {
		
			qui reg `var'
			eststo _all
			matrix coeff = e(b)
			local all_mean = coeff[1,1]
			local all_mean = string(`all_mean', "%9.2f")
			
			gen randXtreat = $randvar*$gottreat 
			qui reg `var' $randvar $gottreat randXtreat
			eststo _second
			matrix coeff = e(b)
			local at_mean = coeff[1,4] + coeff[1,2]		
			local nt_mean = coeff[1,4] + coeff[1,1]	
			local at_mean = string(`at_mean', "%9.2f")
			local nt_mean = string(`nt_mean', "%9.2f")			

			
			qui reg $gottreat $randvar
			eststo _first
			matrix coeff = e(b)
			local comp_frac = coeff[1,1]
			local at_frac = coeff[1,2]
			local nt_frac = 1 - (coeff[1,1] + coeff[1,2])
			local at_frac = string(`at_frac', "%9.2f")
			local nt_frac = string(`nt_frac', "%9.2f")	
			local comp_frac = string(`comp_frac', "%9.2f")				
		
			local `var'_ccm = (`all_mean' - `at_mean' * `at_frac' - `nt_mean' * `nt_frac') / (`comp_frac')
			local `var'_ccm = string(``var'_ccm', "%9.2f")	
			
			display "`var'"
			display `all_mean' "  " `at_mean' "  " `nt_mean' "  " ``var'_ccm' "  " `at_frac' "  " `nt_frac' "  "  `comp_frac'

			qui suest _all _first _second, robust
			qui nlcom dif:  (_b[_all_mean:_cons] - (_b[_second_mean: _cons] + _b[_second_mean: $gottreat ] )*_b[_first_mean:_cons]  -  (_b[_second_mean: _cons] + _b[_second_mean: $randvar ] )* (1 - _b[_first_mean:_cons] - _b[_first_mean:$randvar ] )   )    /   (_b[_first_mean:$randvar ]) 	- _b[_all_mean:_cons], post
			local z = _b[dif]/_se[dif]
			local p = 2*(1-normal(abs(`z')))
			
			loc star "" 
				if `p' < 0.01 {
				loc star "***"
				}
				if `p' <0.05 & `p' > 0.01 { 
				loc star "**"
				}
				if `p' < 0.1 & `p' > 0.05 {
				loc star "*"
				}
				
			local dif_b	 	= string(_b[dif], "%9.2f") + "`star'"
			local dif_se 	= "("+string(_se[dif], "%9.2f")+")"
	
			file write outfile " `: var label `var'' & `all_mean' & ``var'_ccm' & `at_mean' & `nt_mean' & `dif_b' \\  "   _n
			file write outfile " & & & & & `dif_se' \\ "  _n
						
			drop randXtreat	
			
		}		
		file write outfile " \hline " _n	
		file write outfile " Fraction of Sample & 1.00 & `comp_frac' & `at_frac' & `nt_frac' &   \\ "   _n			
		file write outfile " \hline " _n	
		file close outfile
		

			


***
*** Pandemic Panel Table (Table A.9)
***

		clear
		use "$datadir/regdata.dta"
		gen controls = early
		keep if main == 1

		forvalues x = 1(1)6 {
			rename non_hpfl_`x'm non_hpfl_`x'
			rename hp_`x'mo hp_`x'
		}
		keep strata treatment_sum hp_3mo non_hpfl_* HMISID rct_dt early late
		reshape long non_hpfl_ hp_ , i(HMISID) j(month)
		
		gen pandemic = (mofd(rct_dt) + month > mofd(date("3/1/2020","MDY")))
		gen pandemicXtreat = pandemic * treatment_sum
		gen lateXtreat = late * treatment_sum

		foreach var in   non_hpfl_  {

			xi: reg `var'  treatment_sum  late lateXtreat ///
				, cluster(HMISID)
				eststo _panel1
				qui sum `var' if treatment_sum == 0 & late == 0
				estadd scalar cmean = r(mean)

			xi: reg `var'  treatment_sum  pandemic pandemicXtreat ///
				, cluster(HMISID)
				eststo _panel2
				qui sum `var' if treatment_sum == 0 & pandemic == 0
				estadd scalar cmean = r(mean)

			xi: reg `var'  treatment_sum   lateXtreat late i.strata ///
				, cluster(HMISID)
				eststo _panel3
				qui sum `var' if treatment_sum == 0 & late == 0
				estadd scalar cmean = r(mean)

			xi: reg `var'  treatment_sum   pandemicXtreat pandemic i.strata ///
				, cluster(HMISID)
				eststo _panel4
				qui sum `var' if treatment_sum == 0 & pandemic == 0
				estadd scalar cmean = r(mean)
			
		}

		
		*main sample regression table
		esttab	_panel1 _panel3 _panel2 _panel4	using "Tables\table_pandemic_panel.tex",  ///
			replace  frag  nogaps starlevel(* .10 ** .05 *** .01) se ///
			varlabels( treatment_sum "Assigned to Treatment"  lateXtreat "Randomized After X Treatment" late "Randomized After March 2020" pandemicXtreat "Observed After X Treatment" pandemic "Observed After March 2020") ///
			indicate ( "Randomization Strata FE = *strata*") ///
			mtitles("Homeless"  "Homeless"  "Homeless"  "Homeless" ) ///
			stats(cmean r2 N, label("Pre-March Control Mean" "\$R^2\$" "N") ) ///
			b(a2)  alignment(S) nocons


	
***
*** External Validity Table  (Table A.11)
***

	clear
	use "$datadir/regdata.dta"

	keep if (startdate>=td(01jul2019) & startdate<=td(01mar2020))
	
	gen instudy = (rct_dt ~= . & eligible_other == 0)
	
	gen round_score = assmtscorez
	replace round_score = 21 if round_score >= 21 
	
	gen treat = treatment_sum
	replace treat = 0 if assmtscorez >= 8 & assmtscorez <= 13 & treat == .
	replace treat = 1 if eligible_other == 1	
	replace treat = 1 if assmtscorez > 13
	replace treat = 0 if assmtscorez < 8
	
	gen score_cent = (assmtscorez - 10.5)
	gen notstudy = (1-instudy)
	
	
	reg non_hp_6m treat assmtscorez i.strata if instudy == 1, robust
			eststo _reg1
			qui sum non_hp_6m if treat == 0 & instudy == 1
			estadd scalar cmean = r(mean)
			qui sum non_hp_6m if treat == 1 & instudy == 1
			estadd scalar tmean = r(mean)
	
	reg non_hp_6m treat assmtscorez  if instudy == 1, robust
			eststo _reg2
			qui sum non_hp_6m if treat == 0 & instudy == 1
			estadd scalar cmean = r(mean)
			qui sum non_hp_6m if treat == 1 & instudy == 1
			estadd scalar tmean = r(mean)
	
	reg non_hp_6m treat assmtscorez  if instudy == 0, robust
			eststo _reg3
			qui sum non_hp_6m if treat == 0 & instudy == 0
			estadd scalar cmean = r(mean)
			qui sum non_hp_6m if treat == 1 & instudy == 0
			estadd scalar tmean = r(mean)
	
	reg non_hp_6m treat assmtscorez  , robust
			eststo _reg4
			qui sum non_hp_6m if treat == 0 
			estadd scalar cmean = r(mean)
			qui sum non_hp_6m if treat == 1
			estadd scalar tmean = r(mean)	
		
	* regression table
		esttab	_reg1 _reg2 _reg3 _reg4	using "Tables\table_context.tex" ///
		, replace  frag  nogaps starlevel(* .10 ** .05 *** .01) se ///
		varlabels( treat "Assigned to Treatment" assmtscorez "Risk Score"  ) ///
		drop(_cons) ///
		indicate ( "Randomization Strata FE = *strata*") ///
		mtitles("RCT" "RCT" "Not RCT" "All") ///
		stats(cmean tmean r2 N, label("Control Mean" "Treatment Mean" "\$R^2\$" "N") ) ///
		b(a2)  alignment(S) nocons
	


***
***	Main Figures (Figures 1, A.2.a, A.3,
***
	clear
	use "$datadir/regdata.dta"
	
	scalar new_add_max = 0.3
	scalar new_add_inc = 0.05
	scalar new_add_notscc_max = 0.3
	scalar new_add_notscc_inc = 0.05
	scalar shelter_max = 0.1
	scalar shelter_inc = 0.02
	scalar non_hp_max = 0.1
	scalar non_hp_inc = 0.02
	scalar hp_max = 0.7
	scalar hp_inc = 0.1
	scalar dh_max = 0.7
	scalar dh_inc = 0.1
	scalar paymt_tot_max = 4000
	scalar paymt_tot_inc = 500
	scalar paymtany_tot_max = 0.7
	scalar paymtany_tot_inc = 0.1

	foreach samp of varlist early late {
		
		clear
		use "$datadir/regdata.dta"
		keep if `samp' == 1


		* Normalizing variable names so as to loop through creating graphs more easily

		forvalues i = 1(1)12 {
			rename hp_`i'mo hp_`i'
		}

		foreach i in shelter non_hp coord perm rrh tran out dh  {
			forvalues j = 1(1)12 {
				rename `i'_`j'm `i'_`j'
			}
		}

		* making a tempfile in order to re-access it after manipulating the data for the charts below

		
		tempfile maindata
		save `maindata'

		
		***Full non-missing sample
		
			* Creating data for graphs

			*foreach i in new_add any_add shelter non_hp hp coord perm rrh tran out dh {
			foreach i in paymt_tot paymtany_tot new_add new_add_notscc shelter non_hp hp dh  {

				keep HMISID treatment `i'_*
				reshape long `i'_, i(HMISID) j(month)
				rename `i'_ `i'
				drop if month > 12

				tempfile `i'_outcomes
				save ``i'_outcomes'
				
				use `maindata', clear

			}


		* Creating Outcome Graphs for Short Outcome Period of Late Enrollees

			foreach var in paymt_tot paymtany_tot new_add new_add_notscc shelter non_hp hp dh {
				
				if "`samp'" ==	"late" {			
					clear
					local filename ``var'_outcomes'
					use `filename'
					keep if month <= 6
					tab month treatment
					cibar `var', over(treatment month) ///
						vce(cluster HMISID) ///
						graphopts(xtitle("Months Since Random Assignment") ///
						ytitle("Cumulative Proportion") ///
						legend(label(1 "Treatment" 2 "Control")) ///
						scheme(LEO) ylabel(0(`=`var'_inc')`=`var'_max'))  
					graph export "Figures\bar_`var'_`samp'.png", as(png) replace				
				}
				
			}

				
		*** Creating Outcome Graphs for Long Outcome Period of Early Enrollees
		
			if "`samp'" ~= "late" {
			
				clear
				use `maindata'
			
				* Creating data for graphs

			foreach i in paymtany_tot paymt_tot new_add new_add_notscc shelter non_hp hp dh {
					keep HMISID treatment `i'_*
					drop if `i'_12 == .
					reshape long `i'_, i(HMISID) j(month)
					rename `i'_ `i'
					drop if month > 12

					tempfile `i'_outcomes
					save ``i'_outcomes'
					
					use `maindata', clear

				}


			foreach var in paymtany_tot paymt_tot new_add new_add_notscc shelter non_hp hp dh {
					clear
					local filename ``var'_outcomes'
					use `filename'
					cibar `var', over(treatment month) ///
						vce(cluster HMISID) ///
						graphopts(xtitle("Months Since Random Assignment") ///
						ytitle("Cumulative Proportion") ///
						legend(label(1 "Treatment" 2 "Control")) ///
						scheme(LEO) ylabel(0(`=`var'_inc')`=`var'_max') )  
					graph export "Figures\bar_consistsamp_`var'_`samp'.png", as(png) replace
					
					
				}

			}
	
	}


***
*** Funding histogram (Figure A.2.b)
***

	clear
	use "$datadir/regdata.dta"
	
	hist paymt_tot_3 ///
		if early == 1  & treatment_sum == 1 & paymt_tot_3 > 0 ///
		, scheme(LEO) start(0) width(1000)	frac
	graph export "Figures\paymt_distribution_3mo.png", replace
	


***
*** Enrollment over time histogram (Figure A.4)
***

	clear 
	import excel "$datadir\LEO PR-VISPDAT 5.3.2021.xlsx", sheet("Programs") firstrow clear
	rename ClientsUniqueIdentifier HMISID 
	duplicates report HMISID EnrollmentsProjectStartDate ProgramsProjectTypeCode ProgramsFullName
	duplicates drop HMISID EnrollmentsProjectStartDate ProgramsProjectTypeCode ProgramsFullName, force 
	count 

	keep if EnrollmentsProjectStartDate <= date("5/1/2021","MDY")

	format EnrollmentsProjectStartDate %tdMonth_dd,_CCYY
	hist EnrollmentsProjectStartDate  ///
		if  ProgramsProjectTypeCode == "Homeless Prevention" & EnrollmentsProjectStartDate >= date("1/1/2018","MDY") ///
		, width(`=365.25/12') scheme(LEO) freq xlabel(,angle(vertical)) xtitle("") ytitle("Number of Enrollments")

	graph export "Figures\enroll_over_time_3mo.png", replace
	
	
	
***
*** External Validity Figure (Figure A.5)
***
			
	clear
	use "$datadir/regdata.dta"

	keep if (first_assessment_dt>=td(01jul2019) & first_assessment_dt<=td(01mar2020))
	
	gen treat = treatment_sum
	replace treat = 0 if assmtscorez >= 8 & assmtscorez <= 13 & treat == .
	replace treat = 0 if assmtscorez < 8
	replace treat = 1 if assmtscorez > 13
	replace treat = 1 if eligible_other == 1
	
	gen round_score = assmtscorez
	replace round_score = 21 if round_score >= 21 
	
	
	preserve
		gen counter = 1
		collapse (sum) counter (mean) non_hp_6m, by(round_score treat)
			
			twoway 	(scatter non_hp_6m round_score [w=counter] if treat == 0 , msymbol(circle_hollow) ) ///
					(scatter non_hp_6m round_score [w=counter] if treat == 1 , msymbol(diamond_hollow) ) ///
					(pci 0 0 0 0) ///
					(lfit non_hp_6m round_score [w=counter] if treat == 0 , msymbol(circle) ) ///
					(lfit non_hp_6m round_score [w=counter] if treat == 1 , msymbol(diamond) ) ///
					,xtitle("Risk Score") ytitle("% Homeless Services - 6 Mos") scheme(LEO)  ///
					legend(rows(2) order(1 4 2 5)   label(1 "Control") label(2 "Treatment"))


			graph export "Figures/hmis_byscore.png", as(png) replace
					
	restore
	
	
	count

	